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I.  INTRODUCTION 


Reliable  and  accurate  position  information  is  important  when  locating  mobile 
wireless  communication  units.  One  major  commercial  factor  behind  the  recent  interest  in 
localization  is  the  Federal  Communications  Commission  (FCC)  regulation  [1]  that 
requires  wireless  communications  systems  to  provide  Enhanced-911  (E-911)  service. 
There  are  many  other  reasons  where  localization  of  wireless  emitters  is  desirable  or 
essential. 

Our  main  focus  is  the  reduction  in  mean  square  error  for  the  time  difference  of 
arrival  estimate,  which  in  turn  allows  improved  localization.  Wavelet  denoising,  based  on 
the  modified  approximate  maximum  likelihood  (MAML),  the  fourth  order  moment,  or 
their  time  varying  adaptations,  provides  improvement. 

A.  POSITION  LOCALIZATION  SYSTEMS 

A  number  of  position  localization  systems  have  evolved  over  the  years.  Some  of 
these  position  localization  systems  are  Global  Positioning  System  (GPS),  Lor  an  C,  and 
Global  Navigation  Satellite  System  (GLONASS). 

1.  Global  Positioning  System  (GPS) 

GPS  is  the  most  populzir  radio  navigation  aid.  GPS  is  also  used  to  relay  position  of 
cellular  phones  to  public  switched  telephone  networks  (PSTN)  or  to  public  safety 
answering  points  (PSAP).  GPS  obtains  precise  timing  using  a  group  of  satellites  that 
transmit  a  spread  spectrum  L-band  signal  centered  at  1575.42  MHz  [23].  Calculation  of 
the  distances  to  three  different  satellites,  allows  triangulation  to  determine  the  position  of 
the  cellular  phone. 
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2. 


Loran  C 


Loran  C  is  a  navigational  tool  that  operates  in  the  low  frequency  (90-110  kHz) 
band  and  uses  a  pulsed  hyperbolic  system  for  triangulation.  It  has  a  repeatable  accuracy 
in  the  19-90  meters  range  and  is  accurate  to  about  100  meters  with  95  percent  confidence 
and  97  percent  availability.  Like  GPS,  its  performance  depends  on  local  calibration  and 
topography.  GPS  has  replaced  Loran  C  in  most  applications  [13]. 

3.  Global  Navigation  Satellite  System  (GLONASS) 

The  Global  Navigation  Satellite  System  (GLONASS)  is  similar  system  to  GPS. 
GLONASS’s  synchronization  time  is  1/3  of  GPS’s,  i.e.,  less  than  a  minute  [13]. 

B.  METHODS  FOR  LOCATING  CELLULAR  PHONES 

Localization  techniques  can  be  grouped  into  two  categories.  These  are,  position 
localization  systems  that  require  a  modification  of  the  existing  handsets,  and  systems  that 
require  modification  at  the  base  stations.  The  second  category  consists  of  angle  of  arrival 
(AOA),  frequency  difference  of  arrival  (FDOA),  time  of  arrival  (TOA),  or  time 
difference  of  arrival  (TDOA)  estimation  of  the  wavefront  at  the  receiving  platforms. 

1.  Angle  of  Arrival  (AOA) 

A  single  platform  may  be  sufficient  for  AOA  localization  of  a  wireless  transmitter 
on  the  surface  of  the  earth.  Once  the  angle  is  determined,  the  position  of  the  emitter  can 
be  obtained  using  the  intersection  of  the  centerline  of  the  antenna  beam,  with  the  surface 
of  the  earth.  The  signals  to  the  antenna  must  be  coming  from  the  Line-Of-Sight  (LOS) 
direction.  There  is  a  considerable  cost  of  installing  antenna  arrays  and  the  method 
requires  a  relatively  complex  AOA  algorithm  [3, 13]. 

2.  Frequency  Difference  of  Arrival  (FDOA) 
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FDOA  measurements  require  at  least  two  receiving  platforms  and  a  relative 
velocity,  sufficiently  large  so  that  the  differential  Doppler  shift  of  the  two  received 
signals  is  larger  than  the  frequency  measurement  error. 

3.  Time  of  Arrival  (TO  A) 

It  may  be  possible  for  the  base  station  to  indirectly  determine  the  time  that  the 
signal  takes  from  the  source  to  the  receiver  on  the  forward  or  the  reverse  link.  The  total 
time  elapsed  from  the  instant  the  command  is  transmitted  to  the  instant  the  mobile 
response  is  detected,  is  composed  of  the  sum  of  the  round  trip  signal  delay  and  any 
processing  and  response  delay  within  the  mobile  unit.  Due  to  the  variations  in  design  and 
manufacture  of  the  handsets  the  time  estimate  is  difficult  to  obtain.  This  method  is 
susceptible  to  timing  errors  in  the  absence  of  LOS. 

4.  Time  Difference  of  Arrival  (TDOA) 

The  TDOA  estimate  is  taken  as  the  delay,  which  maximizes  the  cross-correlation 
function  between  signals  arriving  at  two  base  stations.  This  also  determines  at  which  base 
station  the  signal  arrives  first.  This  information  yields  a  hyperbolic  localization  curve.  We 
can  localize  the  wireless  transmitter  by  solving  two  hyperbolic  curve  equations. 

It  is  necessary  that  the  code  generators  at  each  receiver  be  synchronized  to  have  a 
common  time  base  [2].  This  form  of  localization  is  useful  in  asynchronous  system  since 
time  of  transmission  need  not  be  known.  To  determine  the  location  of  a  transmitter  in  two 
dimensions,  at  least  three  receivers  are  required.  If  a  major  reflector  effects  the  signal 
components  going  to  all  the  receivers,  the  timing  error  may  get  reduced  in  the  time 
difference  operation. 
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II.  GLOBAL  SYSTEM  FOR  MOBILE  (GSM) 


GSM  is  a  second-generation  cellular  system  standard  that  was  developed  to  solve 
the  fragmentation  problems  of  Europe’s  first  cellular  systems.  GSM  is  the  first  cellular 
system  to  specify  digital  modulation,  network  level  architectures  and  services.  Prior  to 
GSM  it  was  not  possible  to  use  a  given  single  subscriber  unit  throughout  Europe.  GSM 
was  originally  developed  to  serve  as  the  Pan-European  cellular  service  and  promised  a 
wide  range  of  network  service  through  the  use  of  the  Integrated  Services  Digital  Network 
(ISDN).  GSM  is  now  the  world’s  most  popular  standard  for  cellular  and  personal 
communications  equipment. 

A.  GSM  SYSTEM  ARCHITECTURE 

The  GSM  system  architecture  consists  of  three  major  interconnected  subsystems 
that  interact  with  themselves  and  the  users  through  network  interfaces.  The  subsystems 
are  the  Base  Station  Subsystem  (BSS),  Network  and  Switching  Subsystem  (NSS),  and  the 
Operation  Support  Subsystem  (OSS).  The  Mobile  Station  (MS)  is  also  a  subsystem,  but  is 
usually  considered  to  be  part  of  the  BSS  for  architectural  purposes. 

The  BSS  provides  and  manages  radio  transmission  between  the  MS’s  and  the 
Mobile  Switching  Center  (MSC).  The  BSS  also  manages  the  radio  interface  between  the 
MS’s  and  all  other  subsystems  of  GSM. 

The  NSS  manages  the  switching  functions  of  the  system  and  allows  the  MSC’s  to 
communicate  with  other  networks  such  as  the  Public  Switched  Telephone  Network 
(PSTN)  and  Integrated  Services  Digital  Network  (ISDN). 
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The  OSS  supports  the  operation  and  the  maintenance  of  GSM  and  allows  system 
engineers  to  monitor,  diagnose,  and  troubleshoot  all  aspects  of  GSM.  This  subsystem 
interacts  with  the  other  GSM  subsystems. 

1.  GSM  Signal  Specifications 

GSM  utilizes  two  25  MHz  bands,  which  have  been  set  aside  for  system  use  in  all 
member  coimtries.  The  890-915  MHz  band  is  used  for  subscriber-to-base  transmission 
(reverse  link),  and  the  935-960  MHz  band  is  used  for  base-to-subscriber  transmission 
(forward  link).  GSM  uses  Frequency  Division  Duplex  (FDD)  and  a  combination  of  Time 
Division  Multiple  Access  (TDMA)  and  Frequency  Hopped  Multiple  Access  (FHMA) 
schemes  to  provide  stations  with  simultaneous  access  to  multiple  users.  The  available 
forward  and  reverse  frequency  bands  are  divided  into  200  KHz  channels.  These  chaimels 
are  identified  by  their  Absolute  Radio  Frequency  Channel  Number  (ARFCN).  The 
ARFCN  denotes  a  forward  and  reverse  channel  pair,  which  is  separated  in  frequency  by 
45  MHz.  Each  channel  is  time  shared  between  as  many  as  eight  subscribers  using 
TDMA. 

The  radio  transmissions  on  both  forward  and  reverse  link  are  made  at  a  channel 
data  rate  of  270.833  Kbps  using  binary  0.3  OMSK  modulation.  The  bandwidth-bit 
duration  product,  BT,  has  a  level  of  0.3.  The  signaling  bit  duration  is  3.692  ps.  User  data 
is  sent  a  maximum  rate  of  24.7  Kbps.  Each  time  slot  (TS)  has  an  equivalent  time 
allocation  allowing  for  156.25  channel  bits.  From  this,  8.25  bits  of  guard  time  and  6  total 
start  and  stop  bits  are  used  to  prevent  overlap  with  adjacent  time  slots.  Each  TS  has  a 
time  duration  of  576.92  ps,  while  a  single  GSM  TDMA  frame  spans  4.615  ms.  The  total 
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number  of  available  channels  within  a  25  MHz  bandwidth  is  125.  Table  2.1  summarizes 


the  signal  specifications  [5]. 


Parameter 

Specifications 

Reverse  Channel  Frequency 

890-915  MHz 

Forward  Channel  Frequency 

935-960  MHz 

ARFCN  Number 

0  to  124  and  975  to  1023 

Tx/Rx  Frequency  Spacing 

45  MHz 

Tx/Rx  Time  Slot  Spacing 

3  Time  slots 

Modulation  Data  Rate 

270.833333  Kbps 

Frame  Period 

4.615  ms 

Users  per  Frame  (Full  Rate) 

8 

Time  Slot  Period 

576.9  ps 

Bit  Period 

3.692  ps 

Modulation 

0.3  GMSK 

ARFCN  Charmel  Spacing 

200  KHz 

Interleaving  (max  delay) 

40  ms 

Voice  Coder  Bit  Rate 

13.4  Kbps. 

Table  2.1:  GSM  signal  specifications. 


2.  Gaussian  Minimum  Shift  Keying  (GMSK) 

GMSK  is  a  binary  modulation  scheme  which  may  be  viewed  as  a  derivative  of 
Minimum  Shift  Keying  (MSK).  In  GMSK,  the  sidelobe  levels  of  the  spectrum  are  reduced 
by  passing  the  modulating  Non-retum  to  zero  (NRZ)  data  waveforms  through  a  pre¬ 
modulation  Gaussian  pulse-shaping  filter  [12].  Baseband  Gaussian  pulse  shaping 
smoothes  the  phase  trajectory  of  the  MSK  signal  and  hence  stabilizes  the  instantaneous 
firequency  variations  over  time.  This  reduces  the  sidelobe  levels  of  the  transmitted 
spectrum. 

3.  GSM  Channel  Types 

There  are  two  types  of  GSM  logical  channels,  called  traffic  channels  (TCH)  and 
control  channels  (CCH).  Traffic  channels  carry  digitally  encoded  user  speech  or  data  and 
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have  identical  functions  and  formats  on  the  forward  and  reverse  link.  Control  channels 


carry  signaling  and  synchronizing  commands  between  the  base  and  mobile  station. 

4.  Frame  Structure  for  GSM 

There  are  eight  time  slots  (TS)  per  GSM  frame.  The  frame  period  is  4.615  ms.  A 
TS  consists  of  148  bits  which  are  transmitted  at  rate  of  270.833  Kbps  (an  unused  guard 
time  of  8.25  bit  period  is  provided  at  the  end  of  each  burst).  Out  of  the  total  148  bits  per 
TS,  1 14  are  information  bits,  which  are  transmitted  as  two  57  bit  sequences  close  to  the 
beginning  and  end  of  the  burst.  A  26  bit  training  sequence  allows  the  adaptive  equalizer 
in  the  mobile  or  base  station  receiver  to  analyze  the  radio  chaimel  characteristics  before 
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Figure  2.1:  GSM  frame  stmcture. 

decoding  the  user  data.  Stealing  flags  on  the  both  side  of  the  training  sequence  are  used  to 
distinguish  whether  the  TS  contains  voice  or  control  data. 

B.  TRANSMITTER/RECEIVER  SYSTEM 

In  this  section  we  will  look  at  the  GSM  transmitter  and  receiver  and  simulate  some  GSM 
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signals  to  be  used  for  simulation  purposes.  The  overall  structure  of  the  GSM 
transmitter/receiver  system  is  illustrated  in  Figure  2.2. 


Transmitter 


Receiver 
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Channel 

decoder 

/de-interleaver 


De- 

MUX 

De¬ 

modulator 

RF-Rx 

Figure  2.2:Block  diagram  for  the  GSM  transmitter/receiver  system. 


1.  Transmitter  Structure 

The  overall  structure  of  the  transmitter  is  illustrated  in  Figure  2.3.  The  transmitter 
consists  of  four  distinct  functional  blocks. 


Figure  2.3:  The  overall  structure  of  the  transmitter. 
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To  simulate  an  input  data  stream  to  the  channel  encoder/interleaver,  a  sequence  of 
random  data  bits  is  generated.  This  sequence  is  accepted  by  the  MUX,  which  splits  the 
incoming  sequence  to  form  a  GSM  normal  burst.  The  burst  type  requires  that  a  training 
sequence  is  supplied  and  included.  Upon  having  generated  the  prescribed  GSM  normal 
burst  data  structure,  the  MUX  sends  this  to  the  GMSK-modulator.  The  GMSK-modulator 
block  performs  a  differential  encoding  of  the  incoming  burst  to  form  a  NRZ  sequence. 
This  modified  sequence  is  then  subject  to  the  actual  GMSK-modulation  after  which,  the 
resulting  signal  is  represented  as  a  complex  /  and  Q  baseband  signal. 

2.  Receiver  Structure 

The  general  structure  of  the  receiver  is  illustrated  in  Figpre  2.4.  The  demodulator 
accepts  the  GSM  burst,  r,  using  a  complex  baseband  representation.  Based  on  this  data 
sequence,  knowing  the  oversampling  rate  OSR,  the  training  sequence  TRAINING  ,  and 
the  desired  length  of  the  receiving  filter,  Lh,  the  demodulator  determines  a  bit  sequence. 
This  demodulated  sequence  is  then  used  as  the  input  to  the  demultiplexer  (DeMUX) 
where  the  actual  data  bits  are  obtained.  The  remaining  control  bits  and  the  training 
sequence  are  stripped  off  Finally  channel  decoding  and  de-interleaving  is  performed. 


Figure  2.4:  The  overall  structure  of  the  receiver. 
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III.  WAVELETS 


Transforms  are  used  to  obtain  information  that  is  not  observable  in  the  original 
signal.  The  wavelet  transform  provides  a  time-scale  representation  [6].  There  are  other 
transformation,  which  can  give  similar  information,  such  as  the  short  time  Fourier 
transform,  the  Wigner-Ville  distribution,  etc.  This  chapter  will  initially  discuss  Fourier 
and  then  present  wavelet  analysis. 

A.  FOURIER  ANALYSIS 

Fourier  analysis  breaks  a  signal  into  its  constituent  sinusoidal  components.  One 
can  think  of  the  Fourier  analysis  as  a  mathematical  technique  for  transforming  the  signal 
from  a  time-based  to  a  frequency-based  representation  [7]. 

1.  Fourier  Series 

Let  denote  a  periodic  signal  with  period  T^.  By  using  the  Fourier  series 
expansion,  we  can  express  the  signal  as  a  weighted  sum  of  complex  exponentials  [8]: 

,  (3.1) 


where  fo  =  1  /  Tg  and  c„  is  given  by 


(3.2) 


2.  Fourier  Transform 

The  Fourier  transform  of  a  general  continuous  function  g(t)  is  defined  as  [8]: 

G(/)  =  r  .  (3.3) 

J-00 
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G(0  is  a  continuous  function  of  the  variable  f.  The  signal  g(t)  can  be  recovered  by  [8]: 

g(t)  =  r  G{jy^i>df .  (3.4) 

•HOO 

The  functions  g(t)  and  G(f)  are  known  as  a  Fourier  transform  pair. 

Fourier  analysis  is  extremely  useful  because  the  signal’s  frequency  content  is  of 
great  importance,  but  in  transforming  to  the  frequency  domain,  time  information  is  lost. 
However,  many  signals  contain  non-stationary  or  transitory  characteristics  such  as 
trends,  or  abrupt  changes.  These  characteristics  can  be  the  most  important  part  of  the 
signal. 

3.  Short-Time  Fourier  Analysis 

The  short-time  Fourier  transforms  (STFT)  is  an  extension  of  the  Fourier  transform 
designed  to  map  the  signal  into  the  time-frequency  plane.  The  STFT  uses  a  sliding 
window  function  w(t)  to  segment  the  data.  The  STFT  is  given  by 

S{tJ)=  r  g{t)w  ,  (3.5) 

v-oo 

where  w*  {t  - 1)  is  the  sliding  window,  *  represents  conjugation.  The  window  function, 
w(t),  affects  the  characteristics  of  STFT.  As  a  result  of  the  uncertainty  principle,  the  time 
resolution  (Af  )>  and  the  frequency  resolution  (A/’)  of  a  given  signal  are  inversely  related. 
Their  product  has  lower  bound  of  l/4;r ,  which  is  achieved  by  the  Gaussian  window  [9]. 
This  produces  a  trade-off  in  resolution  of  time  and  frequency.  Since  the  choice  of 
window  will  fix  (Ar)and  (A/*)  for  the  entire  time  axis,  the  STFT  partitions  the  time- 
frequency  plane  into  a  uniform  grid.  The  window  can  not  provide  good  time  resolution 
and  good  frequency  resolution  simultaneously  [22]. 

B.  WAVELET  ANALYSIS 
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1. 


Introduction 


A  wavelet  is  an  oscillatory  function  [10].  It  has  its  energy  concentrated  in  time, 
which  allows  the  analysis  of  transient,  non-stationary,  or  time- varying  phenomena. 

2.  The  Continuous  Time  Wavelet  Transform  (CTWT) 

The  continuous  time  wavelet  transform  (CTWT)  is  defined  as  : 


C(T,a)  =  ^  ]g(0 


dt. 


(3.6) 


where  'F(0  is  the  wavelet  function,  ris  the  translation  in  time,  a  denotes  dilation  or 


compression  in  time,  l/^^a  normalizes  the  energy  and  *  denotes  conjugation. 

The  time  and  frequency  resolution  is  controlled  by  the  scale  factor  a.  Low  scales 
correspond  to  high  frequency  wavelets  and  provide  good  time  resolution.  High  scales 
correspond  to  low  frequency  wavelets  with  poor  time  but  good  frequency  resolution.  The 
time-frequency  mapping  of  the  STFT  and  CWT  is  shown  in  Figure  3.1. 

The  STFT  produces  a  uniform  grid  with  a  constant  time  (At)  and  frequency  resolution 
(a/),  while  the  CTWT  has  a  time  and  frequency  resolution  that  depends  on  the  scale. 

3.  The  Discrete  Wavelet  Transform  (DWT) 

Although  the  discretized  continuous  wavelet  transform  enables  the  computation  of  the 
continuous  wavelet  transform  by  computers,  it  is  not  a  true  discrete  transform  [6].  As  a 
matter  of  fact,  the  wavelet  coefficients  are  simply  a  sampled  version  of  the  CWT,  and  the 
information  it  provides  is  highly  redundant,  as  far  as  the  reconstruction  of  the  signal  is 
concerned.  This  redundancy,  on  the  other  hand,  requires  a  significant  amount  of 
computation  time  and  resources. 

The  DWT  provides  sufficient  information,  and  can  offer  a  significant  reduction  in 
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Figure  3.1;  (a)  Frequency-Time  plane  for  SIFT,  (b)  Scale-Time  plane  for  CWT. 

the  computation  time.  It  is  considerably  easier  to  implement  than  the  continuous  wavelet 
transform  and  obtained  by  restricting  the  scale  and  time  parameters  of  the  CWT  to 
discrete  values.  The  DWT  of  a  discrete  signal  g(n)  is  defined  by 

=  1,  (3.7) 

n  \a  \  a  J 

where  a,  b,  and  n  are  the  discrete  versions  of  a,  r ,  and  t  of  Equation  3.6  respectively.  The 
scaling  factor  is  further  restricted  to; 

a  =  al,  J=0,l,...log2(W).  (3.8) 

The  choice  of  UoWill  govern  the  accuracy  of  the  signal  reconstruction  via  the  inverse 
transform.  It  is  popular  to  choose  aQ=2,  since  it  permits  the  implementation  of  fast 

algorithms.  Setting  a  =  2^  produces  octave  bands  called  dyadic  scales.  As  the  scale  level 
is  increased  from  J  to  J+1,  the  analysis  wavelet  is  stretched  in  the  time  domain  by  a  factor 
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of  two.  Hence  the  DWT  output  has  better  frequency  resolution  and  less  precise  time 
resolution  as  the  scale  number  increases. 

If  the  time  shifting  parameter  b  is  restricted  to  kl'^ ,  where  k  is  an  integer,  this 
version  of  the  DWT  is  known  as  the  decimated  DWT  and  can  be  written  as: 

C(2\k2')  =  f^Xg(ri)'¥’{2-^n-k)-  (3.9) 

n=l  v 

where  7  =  0, 1,  2,  •  •  • ,  log^  {N) ,  ^  =  1,  2,  •  •  • ,  NT'^ ,  and  N  is  the  length  of  the  signal  g(n). 

The  term  kl'^  in  the  argument  of  DWT,  indicates  that  C(a,  b)  is  decimated  by  a  factor  of 
two  at  each  successive  scale  J  retaining  only  the  even  points. 

a.  Subband  Coding  and  Multiresolution  Analysis 
The  time-scale  (frequency)  representation  of  the  signal  is  obtained  by 
using  digital  filtering  techniques.  Filters  of  different  cutoff  frequencies  are  used  to 
analyze  the  signal  at  different  scales.  The  signal  is  passed  through  a  series  of  high  pass 
filters  and  low  pass  filters.  Each  one  of  the  filter  is  followed  by  a  two-to-one  decimator. 

The  DWT  analyzes  the  signal  at  different  frequency  bands  with  different 
resolutions  by  decomposing  the  signal  into  a  coarse  and  a  detail  component  at  each  scale. 
The  DWT  employs  two  sets  of  functions,  a  scaling  (low  pass)  and  a  wavelet  (high  pass) 
function.  The  decomposition  of  the  signal  into  different  frequency  bands  is  obtained  by 
successive  highpass  and  lowpass  filtering  of  the  signal  followed  by  the  decimation 
operation. 

This  decomposition  halves  the  time  resolution,  since  the  output  is 
characterized  by  half  the  number  of  samples  compared  to  the  input  signal.  This  operation 
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doubles  the  frequency  resolution,  since  the  frequency  band  of  the  signal  now  spans  only 
half  the  input  frequency  band.  This  procedure,  known  as  subband  coding,  is  repeated, 
b.  WAVELET  ANALYSIS 

Wavelet  transforms  are  very  efficient  and  effective  in  analyzing  a  wide  class  of 
signals  [10]. 

1.  The  wavelet  transform  allows  a  more  accurate  time  description  and 
identification  of  signal  characteristics.  A  wavelet  coefficient  represents  a  component  that 
is  localized  in  time.  The  wavelet  transform  may  allow  a  separation  of  components  of  a 
signal  that  overlap  in  time  or  frequency. 

2.  Wavelets  are  adaptable  and  can  be  designed  to  fit  individual  applications. 
The  discrete  wavelet  transform  is  well  suited  to  digital  implementation. 
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IV.  TIME  DIFFERENCE  OF  ARRIVAL  (TDOA)  ESTIMATION 


TDOA  can  be  employed  to  find  the  position  of  a  GSM  emitter.  Figure  4.1  shows 
a  typical  configuration;  one  emitter  and  a  pair  of  receivers  (discrete  time). 


A.  CORRELATION  FUNCTION 

Frequently  we  would  like  to  know  the  association  between  two  signals,  that  is, 
how  one  signal  is  related  to  the  other.  Correlation  of  signals  is  often  encountered  in  radar 
and  sonar  processing,  and  digital  communications. 

Suppose  that  the  signal  at  receiver  one  is  denoted  by  x(n),  while  y(n)  is  a  time 
shifted  version  of  x(n)  at  receiver  2.  With  additive  noise,  x(n)  zndy(n)  can  be  modeled  as: 

x{n)  =5(«)  +  n,(n)  (4.1) 
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y{n)  =as{n- D)+n2{n)  ,« =  0, 1,  •••,  iV-1  (4.2) 

where  s(n)  is  the  vmknown  source  signal,  «,(«)and  n2(n)are  additive  noise,  D  is  the 
difference  in  arrival  times  at  the  receivers,  a  is  an  attenuation  coefficient,  and  N  is  the 
number  of  samples  in  each  snap  shot  received  at  the  two  receivers. 

The  most  widely  accepted  method  for  obtaining  TDOA  (D  in  Equation  4.2)  uses 
the  cross  correlation  method.  Expectation  of  x(n)  and  y(n)  leads  to 

=  E{a  s(n)  s(n-D-t)+ s(n)n2  (n  -  T)+a  s(n  -  D  -  r)n^  (n) +n^(n)n2(n  -r)} . 

Since  the  noise  and  signal  are  independent,  and  the  noise  has  zero  mean,  the  cross  terms 
of  the  expectation  are  zero.  Also  the  noise  is  independently  distributed,  so  the  expectation 
of  ni(«)  and  «2(n  -  r)  is  also  zero.  The  cross  correlation  of  x(k)  and  y(k)  becomes 

R^(T)=aR,{D  +  r)  .  (4.3) 

Figure  4.2  shows  the  circuit  and  the  block  diagram  for  the  discrete  time  cross 
correlator.  The  cross  correlation  approach  requires  that  receivers  share  a  precise  time 
reference.  The  performance  of  the  TDOA  estimates  can  be  improved  by  increasing  the 
summation  interval.  Once  the  cross  correlation  function  is  computed,  the  value  of  r 
which  maximizes  Equation  4.3  is  used  as  the  estimate  of  the  TDOA.  Figure  4.3  displays 
the  fast  correlation  method.  A  cross  spectral  density  estimate  is  obtained  in  the  fi'equency 
domain,  and  the  cross  correlation  estimate  is  obtained  via  an  inverse  Fourier  transform. 

B.  SIMULATION 

Two  hrmdred  simulations  are  used  to  evaluate  the  performance  of  the  TDOA  method. 

A 

The  z'*  error  is  given  by  e,.  =D-D,-  z  =  l,  2,  •••,A  ,  (4.4) 
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Figure  4.2:  (a)  Detail  of  the  cross  correlator  (b)  Block  diagram  for  the  cross  correlator. 


Figure  4.3:  Fast  cross  correlation  method  using  fast  Fourier  transforms, 
where  N  is  the  number  of  realizations.  The  mean  square  error  (MSB)  is  given  by 

^  1=1 

where  e-  is  the  error  for  the  i”'  realization. 
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V.  WAVELET  DENOISING 


One  important  wavelet  transform  applications  is  noise  reduction.  The  basic  idea 
of  denoising  is  to  retain  the  coefficients  that  preserve  the  signal  while  removing  those  that 
represent  noise.  Two  main  properties  of  the  wavelet  transform  assist  in  separating  the 
noise  coefficients  from  the  rest.  The  first  one  is  that,  by  choosing  the  basis  function  to 
match  the  signal,  the  resulting  decompositions  will  contain  relatively  few  coefficients. 
The  second  one  is  that,  for  a  Gaussian  noise  input,  the  transform  coefficients  will  remain 
Gaussian  [18]. 

Wavelet  shrinkage,  introduced  by  Donoho  and  Johnstone  [18-20],  retains  the 
coefficients  above  a  given  level.  It  discards  the  ones  below  the  given  level.  Denoising 
consists  of  three  steps. 

1 .  T aking  the  wavelet  transform  of  the  input  signal. 

2.  Selecting  a  threshold  and  suppressing  the  noisy  coefficients  by  applying  a 
non-linear  thresholding  technique. 

3.  Performing  the  inverse  wavelet  transform  using  the  modified  coefficients. 
A.  THRESHOLD  VALUES 

The  term  threshold  refers  to  a  constant  that  is  computed.  It  separates  the 
coefficients  that  are  to  be  retained  from  those  that  are  not.  The  noisy  data  can  be  written 
as, 

x{n)=s{n)+CTn{n)  r  =  0,  1,  I,--,  N  ,  (5.1) 

where  s(n)  is  the  input  signal  to  be  recovered  from  the  noisy  data,  n(n)  is  zero  mean,  unit 
variance  additive  white  Gaussian  noise  (AWGN),  N  is  the  length  of  the  signal,  cr  is  the 


19 


standard  deviation  of  the  AWGN.  Algorithms  computing  the  threshold  value  T  require 
estimation  of  cr .  Five  methods  of  computing  thresholds  are  described  below. 

1.  Stein’s  Unbiased  Risk  Estimator  (SURE)  Threshold  ( T^y  ) 

This  method  of  threshold  calculation  was  proposed  by  Donoho  and  Johnstone 
[20].  The  thresholding  is  adaptive:  A  threshold  level  is  assigned  to  each  dyadic  resolution 
level  by  the  principle  of  minimizing  the  Stein  Unbiased  Estimate  of  Risk  (SURE)  [21]  for 
threshold  estimates.  The  SURE  threshold  is  smoothness  adaptive.  If  the  unknown  signal 
contains  jumps,  the  reconstruction  does  also.  If  the  unknown  signal  has  a  smooth 
segment,  the  reconstruction  is  as  smooth  as  the  mother  wavelet  will  allow.  This  statistical 
procedure  calculates  the  estimated  mean  square  error  (risk)  for  a  range  of  threshold 
values,  and  selects  with  the  resulting  minimum  risk. 

2.  Sqtwolog  Threshold  (Tj^) 

Sqtwolog  threshold  uses  a  fixed  form  threshold  yielding  minimax  performance 
multiplied  by  a  small  factor  proportional  to  log(length(signal))  [7] 

=^|2  \og{length{signal))  .  (5.2) 

3.  Heursure  Threshold  (Tj) 

Heursme  threshold  is  mixture  of  SURE  and  sqtwolog  threshold  methods  [7]. 

4.  Minimaxi  Threshold  (TnO 

Minimaxi  threshold  uses  a  fixed  threshold  chosen  to  yield  mimimax  performance 

for  mean  square  error  against  an  ideal  procedure  [7].  The  minimax  principle  is 

used  in  statistics  in  order  to  design  estimators.  It  is  designed  to  select  estimators 

that  minimize  the  worst  case  error  occurring  in  the  set. 
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5. 


Wo-So-Ching  Threshold 


This  threshold  method  was  proposed  by  Wo,  So,  and  Ching.  It  will  be  discussed 
in  detail  in  Chapter  VI. 

The  threshold  techniques,  1  through  4,  are  part  of  Donoho’s  method.  Details  are 
discussed  in  Chapters  VI  and  Vn.  At  the  low  SNR  levels  Donoho’s  method  did  not 
provide  an  improvement.  Threshold  technique  5  was  used  successfully  and  details  are 
discussed  in  Chapter  VI  and  VII. 

B.  THRESHOLDING  (SHRINKAGE) 

Once  a  threshold  value  is  established,  a  number  of  methods  exist  to  apply  the 
threshold  to  suppress  or  modify  the  coefficients  of  the  decomposition.  We  examined 
three  thresholding  methods. 

1.  Hyperbolic  Thresholding 

Hyperbolic  thresholding  was  proposed  by  Wong  and  will  be  discussed  in  detail  in 
Chapter  VI. 

2.  Hard  Thresholding 

Hard  thresholding  can  be  described  as  the  usual  process  of  setting  to  zero  the 
elements  whose  absolute  values  are  lower  than  the  threshold  [7].  For  notational 
convenience,  we  define  x(n)  and  n(n)  as  vectors. 

Z=[jr(0),  41),  42),  4^'-l)] '  „  3, 

Af=[»(0).  4I),  42), -,  4Al-l)r 

Let  Whe  NxN  a  wavelet  transform  matrix.  In  a  vector  form,  the  transformed  output  C  is 
related  to  the  input  vector  X  by  C=W  X ,  where 

C  =  [cij,i),j  =  -W,-,J;  i  =  0,1,2,--, 2^-']  ,  (5.4) 
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and  J  =  logj  {N) .  The  indices  j  and  i  represent  the  scale  and  the  position  in  each  scale, 


while  c(-l,0)  denotes  the  remaining  low-pass  filtered  coefficients. 
The  non-linear  hard  threshold  is  given  by 


c{ij)  = 


ciij)  \for  \c(i,j)\>T 


otherwise  . 


Hard  thresholding  of  the  transformed  coefficients  retains  the  coefficients  that  exceed  the 
threshold  value,  while  all  others  are  set  to  zero. 

3.  Soft  Thresholding 

Soft  thresholding  is  an  extension  of  hard  thresholding.  It  zero  sets  the  elements 
whose  absolute  value  is  lower  than  the  threshold,  and  shrinks  the  remaining  coefficients 
by  the  threshold  value  [7]. 


The  non-linear  soft  thresholding  is  given  by 


V-  -s  ;/or  |c(i,7)|>r 

lo  ;  otherwise. 


The  advantage  of  this  method  is  that  the  results  are  not  as  sensitive  to  the  precise 
value  of  the  threshold  T,  as  in  the  “keep  or  kill”  strategy  of  the  hard  thresholding.  The 
disadvantage  of  this  method  is  that  the  general  shape  of  the  signal  is  modified. 
Thresholding  method  2  and  3  are  part  of  Donoho’s  method,  which  did  not  lead  to 
improvement.  Thresholding  method  1,  hyperbolic  thresholding,  was  used  exclusive  in 
this  work.  Details  can  be  found  in  Chapter  VI  and  VII. 
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VI.  TDOA  ESTIMATION 


Seven  wavelet  based  denoising  methods  are  presented  and  evaluated.  The 
methods  are  denoising  based  on  Donoho’s  method  [7],  denoising  using  the  Wo-So-Ching 
threshold.  [15],  denoising  using  hyperbolic  shrinkage  [16],  denoising  using  median 
filtering  [14],  a  modified  approximate  maximum-likelihood  delay  estimation  based  in 
part  on  [17],  denoising  based  on  the  fourth  order  moment,  and  a  time  varying  adaptation. 
Figure  6.1  is  a  generic  block  diagram  for  all  seven  methods. 

A.  WAVELET  DENOISING  BASED  ON  DONOHO’S  METHOD 

Wavelet  denoising,  proposed  by  Donoho  [7,  18-20],  is  discussed  in  Chapter  V. 
This  method  fails  at  low  SNR’s  and  some  results  are  given  in  Chapter  VII. 

B.  WAVELET  DENOISING  USING  THE  WO-SO-CHING  THRESHOLD 
Prior  to  cross  correlation,  each  one  of  the  sensor  outputs  is  denoised  according  to 

the  Wo-So-Ching  thresholding  rule  to  increase  the  effective  input  SNR.  Wavelet 
denoising  (WD)  is  applied  to  each  received  signal  to  recover  the  corresponding  source 
waveform.  The  restored  signals  are  cross  correlated.  The  TDOA  estimate  is  given  by  the 
argument  at  which  the  cross  correlation  function  attains  its  maximum  value. 

We  define  the  x(n),  y(n),  s(n),  n^(n),  and  «2(")  sequences  in  vector  form  as 

Z=[x(0),x(l),x(2),--,x(iV-l)]" 

J^=[y(0),y(l),y(2),-,y(iV-l)]" 

S=[si0),  5(1),  5(2),  -,  s(N-1)Y  (6.1) 

V,=[«.(0),«i(l),n,(2),-,«,(iV-l)]" 

A^2=h(0),  n,(2),-,n,(N-l)V. 
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WT 


Wavelet  Denoising 


Correlator 


Figure  6.1 ;  System  block  diagram  for  TDOA  estimate  using  wavelet  denoising. 


The  start  time  “0”  denotes  the  first  data  point  in  a  given  block  of  data.  The  transfonned 
output  C  is  related  to  the  input  vector  A"  by  C=}V  X ,  .  where  W  is  a  NxN  orthonoimal 

wavelet  transform  matrix.  C  can  be  decomposed  into  C=Cg  >  where 

Cs  =[cAj,i),  J  =  /  =  0,1,2, •••,2-'-’]  (6.2) 

Q*=KO'>0,  y  =  -l,0,l,-,/;f  =  0,l,2,-,2-'-’]  ,  (6.3) 

are  the  wavelet  transform  coefficients  of  the  source  signal  vector  S  and  the  noise  vector 
Ni^  (k=l,  2),  respectively.  The  indices  i  and  j  represent  the  scale  and  the  position  in  each 

scale,  respectively,  while  (-1,0)  and  (-1,0)  denote  the  low-pass  filtered  coefficients. 

Note  that,  is  still  a  Gaussian  vector.  The  main  idea  of  the  signal  restoration 
using  wavelet  denoising  is  to  adapt  each  c(y,0  to  make  its  value  close  to  c^(J,i)  so  that 
a  good  approximation  of  s(k)  can  be  obtained  after  taking  the  inverse  wavelet  transform. 
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The  Wo-So-Ching  threshold  is  derived  according  to  the  Neyman-Pearson 
criterion  as  it  is  used  in  hypothesis  testing  [15].  This  criterion  is  stated  as  follows: 

Let  c  be  a  Gaussian  random  variable  with  known  variance  cr^ .  Let  a  test  be  conducted 
with  the  following  hypotheses 

=  (6.4) 

and, 

H,-.E{c}^fi,  (6.5) 

We  will  denote  the  acceptance  of  hypothesis  and  as  D^and  ,  respectively.  The 
type  n  error  P{Dq  I H^)  will  be  minimized  for  a  given  when  the  threshold  A  is 

selected  as  follows 

\c-jUQ\<A  =  42cr.eff~\\-a)  ,  (6.6) 

where 

c  is  the  Gaussian  random  variable  with  variance  cr^ , 
is  the  expectation  of  crnider  hypothesis 

A  is  the  Wo-So-Ching  threshold, 
a  is  the  type  I  error  that  is  a=  F(D,  /H„)  and 

etf(y)=~  \e‘'dt  .  (6.7) 

0 

This  denoising  method  discards  the  individual  elements  in  Q  that  are  too  small  in 
magnitude.  The  wavelet  coefficient  c(j,i)  is  regarded  as  totally  due  to  noise  if 
\cU,i)\<^  =  yf2aerf-\l-a)  .  (6.8) 

As  a  result,  (j,  i)  will  be  modified  and  is  given  by 
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if 

otherwise 


(6.9) 


The  denoising  process  is  applied  to  both  x(n)  and  y(n).  The  restored  signals  of  s(n)  are 
denoted  by  j(n)and  a  s{n-D),  respectively.  The  TDOA  estimate  is  given  by  the 
location  of  peak  of  the  cross  correlation  function  of  modified  xfwj  andyfh^  coefficients. 

C.  WAVELET  DENOISING  USING  HYPERBOLIC  SHRINKAGE 


Denoising  of  signals  cormpted  by  noise  uses  three  steps;  wavelet  transform, 
shrinkage,  and  inverse  wavelet  transform.  Shrinkage  reduces  the  value  of  each  coefficient 
in  magnitude  by  an  amount  related  to  the  threshold  value.  The  amount  of  the  shrinkage  is 
determined  by  the  shrinkage  method.  In  this  section  the  hyperbolic  method  is  used. 

Hyperbolic  shrinkage  [16]  is  defined  as: 


c(ij)= 


sign[c{i,  j)]4c{i,jf 

< 

0 

>. 


\c{ij)\  >  X 
|c(i,y)|  <  X ; 


(6.10) 


where  c(i,j)  is  the  detail  function  of  received  signal,  X  is  the  Wo-So-Ching  threshold 
(see  Equation  6.8),  and  c{i,j)  is  the  modified  detail  function  after  shrinkage.  The 
coefficients  are  set  to  the  square  root  of  the  difference  of  the  squares  of  the  values  as  long 
as  the  absolute  value  of  the  coefficient  is  greater  than  the  threshold.  If  the  absolute  value 
of  the  coefficient  is  less  than  the  threshold  then  the  coefficient  is  set  to  zero.  Hyperbolic 
shrinkage  can  give  good  performance  for  a  wide  variety  of  signals.  Other  shrinkage 
techniques  may  offer  better  performance  for  specific  input  signals  and  noise  conditions. 
D.  WAVELET  DENOISING  USING  THE  MEDIAN  FILTER 

In  this  method,  a  median  filter  (filter  length  3)  is  applied  to  the  first  detail  function 
of  the  wavelet  transforms.  The  signals  are  reconstructed  using  the  modified  wavelet 


26 


coefficients. 


The  median  filter  is  a  non-linear  one  dimensional  filtering  technique  that  applies  a 
sliding  window  to  a  sequence  [14].  The  median  filter  replaces  the  center  point  of  the 
window  with  the  median  value  of  all  the  points  contained  in  the  window.  The  length  of 
the  window  is  very  important.  For  example,  for  a  stationary  signal  such  as  a  sinusoid,  a 
longer  window  length  is  better.  But  if  the  signal  is  a  non-stationary,  a  short  window 
length  is  better.  Using  the  median  filter  may  have  drawbacks  if  one  does  not  have  a  priori 
information  about  the  source  signal. 

E.  MODIFIED  APPROXIMATE  MAXIMUM-LIKELIHOOD  DELAY 
ESTIMATION 

An  approximate  maximum  likelihood  (AML)  algorithm  was  proposed  by  Y.  T. 
Chan,  H.  C.  So,  and  P.  C.  Ching  to  estimate  the  TDOA  of  signals  [17].  The  general  idea 
of  this  method  is  to  attenuate  the  frequency  bands  where  the  noise  is  strong  and  to 
enhance  the  frequency  bands  where  the  signal  is  strong.  To  do  this,  we  can  use  the  pre¬ 
filters  as  shown  in  Figure  6.2.  The  pre-filters  i/,(/)  and  chosen  as  follows 


[17] 


Gssif) 


1  ,  GM)  ,  g,(/) 


(6.11) 


Here  (/) ,  (/) ,  and  (/)  denote  the  auto-power  spectra  of  s(k),  n,  {k)  and 

n2 (A:) respectively.  This  choice  of  known  as  maximum  likelihood  (ML) 

weighting.  Multiplying  the  denominator  and  nominator  in  Equation  6.11  by 
G„^„^  if)G„^„^  if)  and  integrate  the  power  spectral  densities,  we  obtain  the  formula  for  the 
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weight  for  each  subband.  This  tends  to  enhance  frequency  bands  where  the  signal  is 
strong. 

The  weights  ,used  as  indicated  in  Figure  6.3,  are  given  by 


*2 


cr 


= 


sd,- 


^n^di^n2dj  '^^n2d,  ) 


^2  /  ^2 


^2 


,i  -1,2,- ■•,J 


(6.12) 


Since  the  noise  density  is  flat  and  the  filter  gain  of  the  wavelet  transform  is  two, 
the  noise  power  will  essentially  remain  same  after  passing  through  each  high  and  low 
pass  filter.  The  noise  power  at  each  subband  can  be  estimated  using  the  formulas  below 


(6.13) 

(6.14) 

=srgmax  R^(t) 

(6.15) 

Figure  6.2:  A  generalized  cross  correlator. 


^2  _  'll 


(6.16) 
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(6.17) 


/v2  ■  ''2 


The  variance  ,  crl^  are  the  noise  powers  at  each  receiver,  while  is  the  signal  power 
at  receiver,  represent  the  noise  powers  at  the  z'*  subband  after  the  wavelet 

transform.  The  signal  power  at  each  subband  can  be  estimated  by 

>  '  =  1.2, --.r  ;  (6,J8) 

k=0  J 

where  N  is  the  length  of  detail  fimction  ,  and  “z”  denotes  the  scale. 

The  denoising  procedure,  shown  in  Figure  6.3,  consists  of  the  three  steps.  The 
modified  subband  sequences  are  obtained  by  weighting  each  subband  by 


d-  =  di  ;  i  =  l,2,---,y 


(6.19) 


The  modified  AML  method  is  applied  to  both  channels.  The  modified  subband 
sequences  d^,dl,---,dl  and  al  are  then  combined  to  form  the  denoised  signals. 

F.  WAVELET  DENOISING  BASED  ON  THE  FOURTH  ORDER  MOMENT 
1.  White  Noise 

The  correlation  fimction  for  white  noise  is  an  impulse  function.  The  samples  of 
white  noise  process  are  uncorrelated.  The  correlation  fimction  of  white  noise  process 
u[n\  has  the  form 


-«o]  • 


(6.20) 


A  white  noise  stationary  process  has  a  correlation  fimction  given  by 


(6.21) 
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SAe^^)=<rl  .  (6.22) 

If  M  is  a  random  vector  of  N  consecutive  samples  from  white  Gaussian  process, 


the  probability  density  function  is  given  by 


1  -ii  A.-1 


=n 


•  i  K — T 

n—0  ^  J,7t^ Q 


2<rl 


(6.23) 


The  samples  of  a  Gaussian  white  noise  process  are  uncorrelated  and  hence  independent. 


2.  Moments  of  the  White  Noise  Process 
The  definitions  of  moments  for  a  stationary  random  process  are; 
=£’{M[n]}=w„  (6.24) 
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M<”U,>E^lnHn+ljhKl!,] 

(6.25) 

m^^^[4,4]=^{m[«]m[« +4]“[^ +4]} 

(6.26) 

(6.27) 

The  first  two  moments  are  the  mean  and  correlation  function  [11].  The  fourth  moment  is 
given  by 

[h.h^h  ]  =  +  A  ]u[n  +  4  ]u[n  + 1,  ]} 

=E^[n]u[n  +  /,  ]}E^[n  +  4  ]w[n  +  4  ]} 

+E\u[n]u[n  +  4  +  /,  ]«[«  +  4  ]} 

+£■{«[« ]«[«  +  4 +  4  ]“[”  +  4  ]}  • 

If =4  =4  =0 

M<^>[0,0,0]  =  3a: 

3.  The  Wavelet  Transform  of  White  Noise 

Since  the  wavelet  transform  is  a  linear  operation,  the  Gaussian  process  remains 
Gaussian  after  passing  through  low  and  high  pass  filters  and  hence  preserves  the  high 
order  moments  properties. 

4.  Fourth  Order  Moment  of  the  Received  Signal 

The  received  signals  are  of  the  form 

x(n)  =  5'(«)+«j(«)  (6.30) 

>»(«)  =  a  ^(/j -Z))+ «2(«)  ,  (6.31) 

where  n,(n)and  n2(n)are  statistically  independent,  zero  mean  Gaussian  random 
variables.  The  fourth  order  moment  estimation  of  the  received  signal  is  given  by. 


(6.28) 


(6.29) 
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(6.32) 


MW[0,0.0]=£{x^[«]) 

=^{[4«]+«i[«]]i 

=i:{s"[n]}+  4£{s^[«]n,[n] }+  6E^\n]r,^[n] } 
+  4i:{s[n]«'[«]  )+£{«,“[«]) 

>E{T^[n]}  =  ^c7l  . 


5.  Mean  and  Standard  Deviation  of  the  Fourth  Order  Moment 


Let  z  =  ,  then  the  mean  of  the  z  becomes 


1  A--!  ,  .  1  N-l 

1=0  1=0 


The  second  moment  of  z  is  given  by. 


.  ,  f  1  A-l  N-\  1 

[•A'  /=0  y=0  J 


1  «  ,  ,) 
N^tS  A=  '  ■  " 


N  N  N'  ’ 


The  variance  of  the  z  is  given  by, 
o]  =E^^}-E^{2} 
=^(9N  +  96)-9<r: 


(6.33) 


(6.34) 


(6.35) 


If  N  is  large,  the  variance  of  z  is  close  to  zero.  We  accoimted  for  the  variance  in  choosing 
a  threshold,  see  Equation  6.42.  Experimentally,  this  amounted  to  a  change  from  to 

3.1c7'  . 

"i 


32 


6.  Denoising 

The  fourth  order  moment  of  a  detail  function,  which  contains  signal  components. 


should  be  greater  than  3cr^, ,  where  ah,  denotes  the  noise  power  at  subband  .  Using 

this  property  the  wavelet  coefficients  which  represent  noise  are  eliminated  and  the  ones 
which  represent  the  signal  are  retained. 

The  noise  power  is  estimated  using  the  formulas  below 


a]  =  arg  max  (r)  (6.36) 

r 

(6.37) 

•  (6.38) 

The  noise  power  after  passing  through  the  first  high  pass  filter  is  given  by 
n  J  =  l,2  ,  (6.39) 


where  G  is  the  gain  of  the  high  and  low  pass  filters.  If  we  generalize  this  formula,  the 
noise  power  of  any  detail  function  t/,.  is 
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After  modifying  all  detail  fimctions,  the  inverse  wavelet  transform  provides  the  denoised 
signal. 

F.  TIME  VARYING  TECHNIQUE 

The  modified  AML  delay  estimation  method  is  based  on  the  idea  of  attenuating 
the  spectral  components  where  the  signal  is  weak.  This  method  works  well  if  the  subband 
information  is  noise  or  signal  only  for  the  duration  of  the  segment. 

Since  some  signals  are  time  varying  or  transient,  the  subbands  can  contain  signal 
only  for  only  part  of  the  segment.  The  method  is  modified  allowing  for  time  variations. 
Each  detail  output  is  divided  into  two  time  blocks.  The  equations  of  section  E  are  used  on 
each  segment. 
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VII.  SIGNAL  DESCRIPTION  AND  SIMULATION  RESULTS 


A.  SIGNAL  DESCRIPTION 

Two  types  of  test  signals  were  used;  a  set  of  generic  signals  and  a  base  band 
GSM  signal.  The  signals  are  described  in  this  section. 

1.  Generic  Signals 

The  13  bit  Barker  code  sequence  was  used  as  the  message  to  create  four  test 
signals.  The  signal  representing  one  code  bit  is  of  the  form  sin(2 *n*  f*nlN)  where 
N-32,  n=0,  1,  ....  21.  Signals  A,  B,  C,  and  D  use  a  value  of  /  of  1,  4,  8,  and  72, 
respectively.  The  generic  signals  are  plotted  in  Figme  7.1. 


(a)  Test  Signal  A 


1 

0.5 
0 

-0.5 
-1 

0  200  400  600 

(d)  Test  Signal  D 


(b)  Test  Signal  B 


0  200  400  600  0  200  400  600 

Figure  7.1:  Generic  signal  (a)  Signal  A,  (b)  Signal  B,  (c)  Signal  C,  (d)  Signal  D. 
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2.  GSM  Signal 

According  to  the  GSM  signal  properties,  see  Chapter  II  ,  a  base  band  GSM 
signal  was  simulated.  The  /  and  0  channel  of  the  GSM  base  band  signal  are  shown  in 
Figure  7.2.  The  MATLAB  code  for  the  GSM  signal  generation  was  adopted  from  [4]. 


(a)  I  channel 


Figure  7.2:  (a)  I  channel  of  GSM  signal  (b)  Q  channel  of  GSM  signal. 

B.  SIMULATION  RESULTS 

The  mean  square  error  (MSB),  as  defined  in  Equation  4.5,  is  used  as  the  criteria  of 
goodness  (i.e.  low  MSB  denotes  a  small  error  and  hence  good  localization).  To  quantify 
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the  global  improvement,  we  compute  the  sum  of  the  MSB’s  (i.e.,  total  error)  for  each 
method  and  compare  it  to  the  total  error  for  the  non-denoised  TDOA  estimate.  The 
improvement  is  given  in  per  cent.  The  total  error  is  obtained  by  summing  the  MSB  over 
the  SNR’s  of  interest  (i.e.,  -3  to  +5  dB). 

1.  Donoho’s  Method 

As  explained  in  Chapter  V,  Donoho  proposed  different  types  of  thresholds  and 
thresholding  techniques.  His  methods  were  evaluated  on  the  generic  signals.  The  mean 
square  error  using  wavelet  denoising  did  not  improve  relative  to  the  non-denoised  version 
for  the  SNR  values  of  interest.  The  mean  square  error  obtained  from  the  correlation  of 
the  raw  signals  is  much  smaller  than  when  using  Donoho’s  method.  Hence  Donoho’s 
method  is  not  investigated  any  further. 

2.  Wavelet  Denoising  Using  the  Wo-So-Ching  Threshold 

a.  Simulation  Results  for  Generic  Signals 

Figure  7.3  (a)  shows  that  there  is  an  slight  improvement  in  the  mean 
square  error  using  Wo-So-Ching  threshold  method  for  generic  signal  A.  As  the  carrier 
frequency  is  increased  (Figure  7.3  (b)  and  Figure  7.4),  the  mean  square  error  also 
increases.  As  long  as  the  number  of  samples  per  binary  bit  is  high  there  is  an 
improvement  using  this  denoising  method,  but  the  advantage  disappears  as  the  carrier 
frequency  increases.  The  results  are  summarized  in  Table  7.1 . 

b.  Simulation  Results  for  the  GSM  Signals 

Figure  7.5  shows  that  there  is  about  a  28%  improvement  in  the  total  mean 
square  error  using  Wo-So-Ching  thresholding  method  for  the  GSM  signal. 
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Figure  7.3:  (a)  MSB  versus  SNR  with  and  without  denoising  using  the  Wo-So-Ching 


thresholding  method  for  generic  signal  A.  (b)  MSB  for  generic  signal  B. 


(a) 


Figure  7.4:  (a)  MSB  versus  SNR  with  and  without  denoising  using  the  Wo-So-Ching 


threshold  method  for  generic  signal  C.  (b)  MSB  for  generic  signal  D. 
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SNR 

Xcorr 

(A) 

Wo-So-Ching 

threshold 

(A) 

Xcorr 

(B) 

Wo-So-Ching 

threshold 

(B) 

5 

0.385 

0.25 

0 

0 

4 

0.415 

0.275 

0 

0 

3 

0.485 

0.425 

0 

0 

2 

0.715 

0.535 

0 

0 

1 

0.795 

0.78 

0.32 

0.32 

0 

0.805 

0.725 

0 

0.96 

-1 

1.145 

1.02 

2.89 

3.84 

-2 

1.62 

1.34 

4.505 

6.49 

-3 

2.085 

1.815 

8.38 

12.36 

Table  7.1:  (a)  MSE  versus  SNR  for  with  and  without  denoising  using  the  Wo-So-Ching 
threshold  method  for  generic  signal  A  and  B. 


SNR 

Xcorr 

(C) 

Wo-So-Ching 
threshold 
_ CC) 

Xcorr 

(D) 

Wo-So-Ching 

threshold 

(D) 

5 

0 

0 

0 

0 

4 

0 

0.16 

0 

0 

3 

0 

0.56 

0 

0 

2 

0.08 

0.64 

0 

0 

1 

1.2 

3.52 

0 

0 

0 

2.72 

5.84 

0.32 

1.92 

-1 

3.52 

8.56 

1.6 

5.12 

-2 

5.92 

16.56 

5.095 

9.7 

-3 

8.88 

28.08 

5.665 

15.62 

Table  7.1:  (b)  MSE  versus  SNR  for  with  and  without  denoising  using  the  Wo-So-Ching 
threshold  method  for  generic  signal  C  and  D. 
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Figure  7.5;  MSB  versus  SNR  with  and  without  denoising  using  the  Wo-So-Ching 
threshold  method  for  the  GSM  signal. 

3.  Wavelet  Denoising  Using  Hyperbolic  Shrinkage 
a.  Simulation  Results  for  Generic  Signals 

Figure  7.6  shows  that  there  is  an  improvement  in  the  mean  square  error 
using  hyperbolic  shrinkage  method  for  generic  signal  A  and  B.  Again  as  the  carrier 
frequency  is  increased  (Figure  7.8),  the  mean  square  error  also  increases.  The  results  are 
summarized  in  Table  7.2.  Comparing  Table  7.1  and  Table  7.2,  we  can  conclude  that  the 
hyperbolic  shrinkage  method  provides  better  result  than  the  Wo-So-Ching  threshold 
method. 
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b. 


Simulation  Results  for  the  GSM  Signal 


Figure  7.8  shows  that  there  is  about  a  48  %  improvement  in  the  total 
mean  square  error  using  Hyperbolic  Shrinkage  method  for  GSM  signal.  Comparing 
Figure  7.5  and  7.8  we  can  conclude  that  hyperbolic  shrinkage  method  provides  better 
result  for  the  GSM  signal  than  Wo-So-Ching  thresholding  method. 


(a) 


Figure  7.6:  (a)  MSB  versus  SNR  with  and  without  denoising  using  the  Hyperbolic 
Shrinkage  method  for  generic  signal  A.  (b)  MSB  for  generic  signal  B. 
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SNR 


Figure  7.7:  (a)  MSB  versus  SNR  with  and  without  denoising  using  the  Hyperbolic 


Shrinkage  method  for  generic  signal  C.  (b)  MSB  for  generic  signal  D. 


SNR 

Xcorr 

(A) 

Hyperbolic 

shrinkage 

(A) 

Xcorr 

(B) 

Hyperbolic 

shrinkage 

(B) 

0.375 

0.2000 

0 

0 

0.53 

0.3200 

0 

0.08 

0.52 

0.4100 

0 

0.08 

0.625 

0.5400 

0 

0.16 

1 

0.96 

0.7350 

0.32 

0.24 

0 

1.07 

0.8200 

0 

1.12 

-1 

1.345 

1.0050 

2.89 

2.8 

1.36 

1.2100 

4.505 

3.6 

H 

1.895 

1.3300 

8.38 

7.84 

■9 

2.425 

2.2400 

12.025 

10.465 

mm 

2.57 

3.5350 

38.24 

17.185 

191 

87.41 

14.5850 

268.01 

190.41 

Table  7.2:  (a)  MSB  versus  SNR  for  with  and  without  denoising  using  Hyperbolic 
shrinkage  method  for  generic  signals  A  and  B. 
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SNR 

Xcorr 

(C) 

Hyperbolic 

shrinkage 

(C) 

Xcorr 

(D) 

Hyperbolic 

shrinkage 

(D) 

5 

0.1200 

0 

0 

4 

0.3200 

0 

0.24 

3 

0.4800 

0 

0.32 

2 

0.16 

1.0000 

0 

1.76 

1 

1.04 

2.0400 

0 

2.8 

0 

1.6 

3.2200 

0 

4.015 

-1 

4.08 

4.3800 

1.645 

7.625 

-2 

5.0800 

0.96 

12.735 

-3 

9.76 

10.6450 

8.76 

87.35 

-4 

14.8 

16.6750 

15.385 

29.055 

-5 

24.96 

20.1500 

22.97 

680.8 

-6 

118.8 

163.7250 

57.635 

1278.2 

Table  7.2:  (b)  MSE  versus  SNR  for  with  and  without  denoising  using  Hyperbolic 
shrinkage  method  for  generic  signals  C  and  D. 


Figure  7.8:  MSE  versus  SNR  with  and  without  denoising  using  the  Hyperbolic  Shrinkage 
method  for  the  GSM  signal. 
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4. 


Wavelet  Denoising  Using  the  Median  Filter 


a.  Simulation  Results  for  Generic  Signals 

Figiire  7.9  shows  that  there  is  an  improvement  in  the  mean  square  error 
using  median  filtering  on  generic  signals  A  and  B.  As  the  carrier  frequency  is  increased 
(Figure  7.10),  the  mean  square  error  also  increases.  The  results  are  summarized  in  Table 
7.3.  Comparing  Table  7.1,  7.2  and  7.3  we  can  conclude  that  median  filtering  has  the  best 
result  for  the  generic  signal  B  but  for  the  other  signals  the  Hyperbolic  shrinkage  method 
provides  better  results. 


(a) 


Figxure  7.9:  (a)  MSB  versus  SNR  with  and  without  denoising  using  the  median  filtering 
method  for  generic  signal  A.  (b)  MSB  for  generic  signal  B. 
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(a) 


(b) 


Figure  7.10:  (a)  MSE  versus  SNR  with  and  without  denoising  using  the  median  filtering 
method  for  generic  signal  C.  (b)  MSE  for  generic  signal  D. 
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Table  7.3:  (a)  MSE  versus  SNR  for  with  and  without  denoising  using  the  median  filtering 
method  for  generic  signals  A  and  B. 
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Table  7.3:  (b)  MSB  versus  SNR  for  with  and  without  denoising  using  the  median  filtering 
method  for  generic  signals  C  and  D. 

b.  Simulation  Results  for  the  GSM  Signal 
Figure  7.11  shows  that  there  is  about  a  42%  improvement  in  the  total 
mean  square  error  using  the  median  filtering  method  for  the  GSM  signal.  Comparing 
Figure  7.5,  7.8  and  7.11  we  can  conclude  that  hyperbolic  shrinkage  method  gives  the  best 
result  for  the  GSM  signal. 

5.  Modified  Approximate  Maximum-Likelihood  Delay  Estimation 
a.  Simulation  Results  for  Generic  Test  Signals 

Figure  7.12  and  7.13  show  that  there  is  a  significant  improvement  in  the 
mean  square  error  for  all  generic  signals  when  using  the  modified  AML  estimation 
method.  The  results  are  summarized  in  Table  7.4.  Comparing  Table  7.1,  7.2,  7.3  and  7.4 
we  can  conclude  that  the  modified  AML  estimation  method  gives  the  best  result. 
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Figure  7.11:  MSB  versus  SNR  with  and  without  denoising  using  the  median  filtering 
method  for  the  GSM  signal. 


(a) 


Figtire  7.12:  (a)  MSB  versus  SNR  with  and  without  denoising  using  the  modified  AML 
estimation  method  for  generic  signal  A.  (b)  MSB  for  generic  signal  B. 
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SNR 


Figure  7.13:  (a)  MSB  versus  SNR  with  and  without  denoising  using  the  modified  AML 
estimation  method  for  generic  signal  C.  (b)  MSB  for  generic  signal  D. 
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Table  7.4:  (a)  MSB  versus  SNR  for  with  and  without  denoising  the  modified  AML 
estimation  method  for  generic  signals  A  and  B. 
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Table  7.4:  (b)  MSE  versus  SNR  for  with  and  without  denoising  using  the  modified  AML 
estimation  method  for  generic  signals  C  and  D. 


b.  Simulation  Results  for  the  GSM  Signal 

Figure  7.14  shows  that  there  is  about  a  79%  improvement  in  the  total 
mean  square  error  using  the  modified  AML  estimation  method  relative  to  the  undenoised 
version.  Comparing  Figure  7.5,  7.8,  7.11  and  7.14  we  can  conclude  that  the  modified 
AML  estimation  method  provides  the  best  result  for  the  GSM  signal. 

6.  Wavelet  Denoising  Based  on  the  Fourth  Order  Moment 
a.  Simulation  Results  for  Generic  Signals 

Figure  7.15  and  7.16  show  that  there  is  a  significant  improvement  in  the 
mean  square  error  using  wavelet  denoising  based  on  the  fourth  order  moment  method  for 
all  generic  signals.  The  results  are  summarized  in  Table  7.5.  Comparing  this  method  with 
the  other  methods,  we  can  conclude  that  wavelet  denoising  based  on  the  fourth  order 
moment  is  the  second  best  method. 
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Figure  7.14:  MSB  versus  SNR  with  and  without  denoising  using  the  modified  AML 


estimation  method  for  the  GSM  signal. 
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(a) 


Figure  7.16:  (a)  MSB  versus  SNR  with  and  without  denoising  using  the  wavelet 


denoising  based  on  the  fourth  order  moment  method  method  for  generic  signal  C.  (b) 


MSB  for  generic  signal  D. 
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Table  7.5:  (a)  MSB  versus  SNR  for  with  and  without  denoising  wavelet  denoising  based 
on  the  fourth  order  moment  method  for  generic  signals  A  and  B. 
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Table  7.5:  (b)  MSE  versus  SNR  for  with  and  without  denoising  using  wavelet  denoising 
based  on  the  fourth  order  moment  method  for  generic  signals  C  and  D. 

b.  Simulation  Results  for  the  GSM  Signal 

Figure  7.17  shows  that  there  is  about  a  63%  improvement  in  the  mean 
square  error  using  wavelet  denoising  based  on  the  fourth  order  moment  method. 
Comparing  Figure  7.5,  7.8,  7.11,  7.14  and  7.17  we  can  conclude  that  wavelet  denoising 
based  on  the  fourth  order  moment  ranks  second  in  performance  for  the  GSM  signal. 

7.  Time  Varying  Technique 

As  explained  in  Chapter  VI,  we  can  address  signals  that  are  non-stationary  (i.e., 
the  signal  does  not  stay  in  a  given  frequency  band  for  the  length  of  the  segment  or  has 
modulation  characteristics).  In  principle,  in  any  given  subband  we  try  to  find  several 
weights.  This  allows  us  to  keep  the  information  for  the  frequency  band  when  the  signal  is 
strong.  Figure  7.18  allows  a  comparison  of  the  modified  AML  method  with  the  time 
varying  version  in  which  the  segments  are  partitioned  into  two  parts  for  every  scale.  The 
improvement  in  the  total  mean  square  error  using  the  time  varying  modified  AML 
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method  is  about  81%.  There  is  no  drastic  improvement  relative  to  the  modified  AML,  but 
we  believe  that  if  the  signal  has  a  time  varying  property,  the  new  method  based  on  the 
time  varying  technique  will  give  better  results. 


SNR 


Figure  7.17:  MSB  versus  SNR  with  and  without  denoising  using  the  wavelet  denoising 
based  on  the  fourth  order  moment  method  for  the  GSM  signal. 
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Figure  7.18:  MSE  versus  SNR  with  and  without  denoising  using  the  time  varying 
modified  AML  and  modified  AML  technique  on  the  GSM  signal. 
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VIIL  CONCLUSION 


A.  Summary 

The  accuracy  of  algorithms  to  localize  wireless  communication  units,  was  studied. 
Chapter  I  reviewed  the  need  to  locate  cellular  transmitters  and  some  of  the  existing 
localization  systems.  The  time  difference  of  arrival  (TDOA)  method  was  used  to  locate 
emitters.  Estimation  of  the  TDOA  was  based  on  the  cross  correlation  function.  To 
increase  the  accuracy  of  the  TDOA  estimation,  the  noise  in  the  received  data  was  reduced 
(denoising).  The  wavelet  transform  was  employed  to  minimize  the  noise.  Several 
denoising  methods  were  examined. 

The  MSE  for  all  methods  for  the  GSM  signal  is  plotted  in  Figure  8.1.  The 
modified  AML  delay  estimation  method  provided  the  best  result,  wavelet  denoising  using 
the  fourth  order  moment  method  ranks  second.  All  methods  provide  better  results  than 
the  correlation  of  the  raw  time  domain  signals.  The  time  varying  technique  outperform 
the  modified  AML  method  for  SNR  values  greater  or  equal  to  -2  dB.  We  believe  that  if 
the  signal  has  time  varying  properties,  the  time  varying  adaptation  will  be  the  superior 
method. 

For  the  GSM  signal  the  probability  of  no-error  is  plotted  in  Figure  8.2  and  8.3. 
The  Wo-So-Ching,  the  hyperbolic  shrinkage  and  the  median  filtering  techniques,  see 
Figure  8.2,  provide  worse  results  than  the  direct  correlation  of  the  raw  time  domain 
signals.  In  these  techniques,  a  threshold  value  was  calculated  and  then  according  to  the 
threshold,  each  coefficient  was  modified.  These  three  techniques  can  modify  the 
coefficients,  which  represent  the  signal  in  the  subband.  Hence  we  reduced  the  mean 
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square  error  by  using  these  three  techniques  but  we  also  decreased  the  probability  of  no¬ 
error. 

The  modified  AML,  the  fourth  order  moment  and  the  time  varying  AML 
techniques,  see  Figure  8.3,  provide  better  results  than  the  direct  correlation  of  the  raw 
time  domain  signals.  These  techniques  keep  the  subbands  or  portions  of  it,  which 
represent  the  signal.  The  subband,  which  represents  noise  only,  is  eliminated.  The  time 
varying  AML  method  produces  the  highest  probability  of  no-error. 


Figure  8.1 :  Plot  of  all  denoising  methods  for  the  GSM  signal. 
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Figure  8.2:  Plot  of  the  probability  of  no-error  for  the  Wo-So-Ching  threshold,  Hyperbolic 
Shrinkage  and  Median  filtering  techniques. 


B.  Recommendations 

1 .  The  work  was  done  exclusively  using  a  base  band  GSM  signal.  A  follow 
on  study  should  address  GSM  signals  at  the  IF  or  RF  level,  to  assess  potential 
improvement  of  the  MSB  and  hence  in  localization. 

2.  White  Gaussian  noise  was  used  in  the  simulations.  The  effects  of  fading 
signals  should  be  included  in  a  follow  on  study. 
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probability  of  no-error 


SNR 


Figure  8.3:  Plot  of  the  probability  of  no-error  for  the  Time  varying  AML,  Fourth  order 


moment  and  Median  filtering  techniques. 
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